use std::mem::swap;

use rand::random;

/// greatest common divisor of a and b
fn gcd(a: u128, b: u128) -> u128 {
    if b == 0 {
        return a;
    }
    gcd(b, a % b)
}
/// |a-b|
fn dis(a: u128, b: u128) -> u128 {
    if a > b {
        a - b
    } else {
        b - a
    }
}

/// calculate (a * b) % mo without overflow
fn qmul(mut a: u128, mut b: u128, mo: u128) -> u128 {
    if a < u64::MAX as u128 && b < u64::MAX as u128 {
        return a * b % mo;
    }
    if a < b {
        swap(&mut a, &mut b);
    }
    let mut res = 0;
    while b > 0 {
        if (b & 1) == 1 {
            res = (res + a) % mo;
        }
        a = (a + a) % mo;
        b >>= 1;
    }
    res
}

/// quickly calculate (a ^ p) % mo
fn qpow(mut a: u128, mut p: u128, mo: u128) -> u128 {
    let mut res = 1;
    while p > 0 {
        if (p & 1) == 1 {
            // ans = ans * a % mo;
            res = qmul(res, a, mo);
        }
        // a = a * a % mo;
        a = qmul(a, a, mo);
        p >>= 1;
    }
    res
}

/// use Miller-Rabin algorithm
fn is_prime(n: u128) -> bool {
    if n < 3 || (n & 1) == 0 {
        return n == 2;
    }
    if n % 3 == 0 {
        return n == 3;
    }
    let (mut u, mut t) = (n - 1, 0);
    while (u & 1) == 0 {
        u /= 2;
        t += 1;
    }
    for _ in 0..10 {
        let a = random::<u128>() % (n - 3) + 2;
        let mut v = qpow(a, u, n);
        if v == 1 {
            continue;
        }
        let mut s = 0;
        for _ in 0..t {
            if v == n - 1 {
                break;
            }
            // v = v * v % n;
            v = qmul(v, v, n);
            s += 1;
        }
        if s == t {
            return false;
        }
    }
    true
}

/// f(x) = x^2 + c  (% mo)
fn f(x: u128, c: u128, mo: u128) -> u128 {
    (qmul(x, x, mo) + c) % mo
}

/// find any divisor of x
fn poletard_rho(x: u128) -> u128 {
    let mut t: u128 = 0;
    let c: u128 = random::<u128>() % (x - 1) + 1;
    let mut s = t;
    let mut goal = 1;
    let mut val = 1;
    loop {
        for step in 1..=goal {
            t = f(t, c, x);
            val = qmul(val, dis(t, s), x);
            if val == 0 {
                return x;
            }
            if step % 127 == 0 {
                let d = gcd(val, x);
                if d > 1 {
                    return d;
                }
            }
        }
        let d = gcd(val, x);
        if d > 1 {
            return d;
        }
        goal <<= 1;
        s = t;
        val = 1;
    }
}

/// deep first search the maximum divisor
fn dfs(mut x: u128, current_ans: u128) -> u128 {
    if x <= current_ans || x == 1 {
        return current_ans;
    }
    if is_prime(x) {
        return x.max(current_ans);
    }
    let d = poletard_rho(x);
    while x % d == 0 {
        x /= d;
    }
    let t = dfs(x, current_ans);
    dfs(d, t.max(current_ans))
}

pub fn find_max_prime_factor(number: u128) -> u128 {
    dfs(number, 1)
}
